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k^ ' Abstract. We consider the evolution in full general relativity of a family of linearly unstable 

O , isolated spherical neutron stars under the effects of very small, perturbations as induced by the 

\^ ^ , truncation error Using a simple ideal-fluid equation of state we find that this system exhibits 

' a type-I critical behaviour, thus confirming the conclusions reached by LiebHng et al. (T) for 

Cn ' rotating magnetized stars. Exploiting the relative simphcity of our system, we are able carry 

^Nj ' out a more in-depth study providing solid evidences of the criticality of this phenomenon and 

also to give a simple interpretation of the putative critical solution as a spheiical solution with 
the unstable mode being the fundamental F-mode. Hence for any choice of the polytropic 
constant, the critical solution will distinguish the set of subcritical models migrating to the 
stable branch of the models of equihbrium from the set of subcritical models collapsing to a 
t' , ' black hole. Finally, we study how the dynamics changes when the numerically perturbation 

{2JQ' is replaced by a finite-size, resolution independent velocity perturbation and show that in such 

cases a nearly-critical solution can be changed into either a sub or supercritical. The work 
reported here also lays the basis for the analysis carried in a companion paper, where the 
^Nj ' critical behaviour in the head-on collision of two neutron stars is instead considered (l)- 

> 
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QQ , PACS numbers: 04.25.Dm, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60.Jd 
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t^^ ■ 1. Introduction 

o. 

^D . Critical phenomena in general relativity were first discovered by Choptuik f3l| for the 

gravitational collapse of a massless scalar field. After his seminal work these phenomena 

were also discovered for a wide range of systems, including massive scalar fields and ultra- 

k> ■ relativistic fluids ll4 HT2l (see also ifTSJI for a recent review). Our interest here is to further 

''T^ . develop the analysis of critical phenomena in the presence of (perfect) fluids in general and 

C^ ' of relativistic spherical stars in particular. 

We recall that the first simulations of critical collapse of a perfect fluid were performed by 
Evans and Coleman 14|. They used an ultra-relativistic equation of state (EOS) with polytropic 
exponent r = 4/3 (radiation) and found evidence of a continuously self-similar (CSS) critical 
solution. Following works by PI and fT4'| showed that type-Il critical phenomena occur for 
every value of F between 1 and 2 with a CSS critical solution. The ultra-relativistic EOS, 
that is one in which the pressure is proportional to the energy density p = (F — l)e, is the 
only scale-invariant equation of state in general relativity, thus the only one compatible with 
the existence of an homothetic vector field. In other words, fluids obeying an ultra-relativistic 
EOS are the ones that admit self-similar solutions in general relativity ifTSJI . Nevertheless 
evidence of CSS critical solutions were found by |]6| also for perfect fluids with an ideal- 
gas EOS. They conjectured that, as type-II phenomena are kinetically dominated, an ideal 
gas would behave like an ultra-relativistic gas and thus nearly-critical solutions will approach 
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the corresponding ultra-relativistic CSS solution. They were actually able to show that this 
is indeed the case by comparing ultra-relativistic CSS solutions, obtained using the self- 
similarity ansatz, to the corresponding ideal-gas solution. This suggested the possibility of 
observing type II critical phenomena in the context of neutron star collapse. 

The first studies in this direction were performed by Novak I?). He studied solutions 
with velocity induced perturbation of a Tolman-Oppenheimer-Volkoff (TOV) background and 
found a type-II critical phenomena using a stiff ideal gas EOS {i.e. with F = 2) and even a 
more realistic tabulated EOS. Noble [JS] studied critical phenomena in neutron star collapse 
in great detail by using the gravitational interaction with a massless scalar field to perturb 
spherical-star configurations. He confirmed the observation by Novak that a minimum mass 
is required to trigger collapse and also found that, for very high mass spherical stars, type-I 
critical phenomena is observed with each model oscillating around a corresponding solution 
on the unstable branch. More recently [J9] performed an accurate study using the same setup 
of Q and was able to show that, even for initial data of spherical stars, the scaling exponent 
in the mass relation law is compatible with the exponent in the corresponding ultra-relativistic 
case. 

Type-I critical phenomena in neutron star collapse have seen a renewed interested mainly 
due to the discover by ifTol of critical phenomena in head-on collision of neutron stars (NS). 
They considered families of two equal mass NS, modeled with an ideal gas EOS, boosted 
towards each other at a prescribed speed and varied the mass of the stars, their separation, 
velocity and F-parameter They found that at the threshold of black hole formation a type- 
I critical phenomena can be observed, with the putative solution being an oscillating star. 
In a successive paper 1161 repeated the same results by performing head-on collision of 
Gaussian packets. They also claimed that the critical solution is a new kind of metastable 
object and not a perturbed spherical star as its mass is less then the maximum allowed mass 
of a spherical star for the same EOS [[T0| and its oscillation frequencies are one order of 
magnitude larger then the eigenfrequencies of a spherical star with the same total baryonic 
mass lfT6l . This claims have partly been rejected by the work carried out in [2|, which also 
investigated critical phenomena in head-on collision of equal-mass spherical stars and found 
that the critical solution can be compared to a spherical star solution sitting on the unstable 
branch. 

Motivated by these recent results we have investigated phase transitions in families of 
linearly unstable spherical stars, a case substantially unexplored in previous systematic works. 
We show that this transition is a critical one and that the putative critical solution can be 
interpreted as an oscillating spherical star We further study the effects of the introduction of 
velocity-induced perturbations, as the ones used by Q, on nearly critical solutions. 

The remainder of this paper is organized as follows. In section |2] we give a brief 
introduction to critical phenomena in general relativity. In section[3]we describe the numerical 
settings of the simulations and the properties of the used initial data. In section |4] we show 
in details our results, while section |5] is dedicated to conclusions and discussion. We use a 
spacetime signature (— , +, +, +), with Greek indices running from to 3 and Latin indices 
from 1 to 3. We also employ the standard convention for the summation over repeated 
indices. Unless otherwise stated, all the quantities are expressed in a system of units in which 
c = G = Mq = 1. 

2. Critical phenomena in gravitational collapse 

In what follows we give a brief overview of critical phenomena in gravitational collapse and 
which will be useful to cast our results in the more general context of critical phenomena in 
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general relativity. We refer the interested reader to |fT3]| for a more systematic presentation. 

2.1. Self-similarity 

Before dwelling on critical phenomena and because self-similarity plays a central role in 
this context, it is useful to recall briefly the definitions of "continuous" self-similarity and 
"discrete" self-similarity. We refer the interested reader to ifTTl for a more detailed discussion. 
We recall that a spacetime is said to be continuously self-similar if there exist a 
vector field, ^^ such that V(^^,y) ~ g^,^. Vector fields satisfying this condition are said 
to be "homothetic" as we can easily construct a one-parameter group of transformations, 
(j)s:x^ I— )• ?/^(s), where j/^(s) is the integral curve associated with ^^ passing through a;^. 
It is then easy to see that (ps is an homothetic transformation as the associated push-forward, 
acts as a rescaling on the metric 

For this reason in a system of coordinates adapted to the self-similarity 

the metric coefficients read 

9tj.u{r,x'') = e"^^5^i,(x*) , (3) 

and the new metric g^y appears explicitly self-similar, i.e. independent of the r. 

Similarly, a spacetime is said to be "discretely self-similar" (DSS) if a discrete version 
of ([T]) holds. In particular, Gundlach ||T8| defines a spacetime to be DSS if there exist a 
diffeomorphism and a real constant A such that for any positive integer n 

/ 1^ \n 2nA ^a\ 

(0 ) 9tiu = e g^u ■ (4) 

In coordinates adapted to the self-similarity a point P with coordinates (r, x*) is mapped by 
cf> into (r — A, a;*) and the metric can be written as 

9tj.i^{T, x') = e^'^^g^y{T, x*) , (5) 

where 

g^,.(T + A,a;') =.gp^(r, a:'). (6) 

Thus, if V'V is timelike and induces a Cauchy foliation of the spacetime, we can give a 
physical interpretation of the dynamics of DSS solutions as a combined effect of a rescaling 
and a periodic "echoing" of the geometry. 

2.2. The basic concepts 

Let us consider a group of one-parameter families of solutions, S[P], of the Einstein equations 
such that for every P > P*, S[P] contains a black hole and for every P < P*, S[P] is a 
solution not containing singularities. We say that these families exhibit a critical phenomenon 
if they have the common property that, as P approaches P*, S[P] approaches a universal 
solution 5[P'^], i.e. not depending on the particular family of initial data, and that all the 
physical quantities of 5[P] depend only on |P — P*|. In analogy with critical phase transitions 
in statistical mechanics, these phenomena are then classified as type-II or type-I critical 
phenomena lfT3l . In what follows we briefly recall the differences between the two classes. 
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P> P* 



Figure 1. Phase space picture of type-II critical phenomena. The surface C represents the 
critical manifold, separating the basins of attraction of A and B. The line 7 represents a 
generic one-parameter family of initial data intersecting the critical manifold in P* . Generic 
initial data starting at Z{0) will evolve towards Aot B following the arrows Z{t), data near 
the threshold will be marginally attracted towards the critical solution Z* . Points exactly on 
the critical manifold will be attracted to the critical solution. 



2.2.1. Type-II critical phenomena. Type-II critical phenomena involve the existence of a 
CSS or DSS solutions sitting at the threshold of black-hole formation. They are characterized 
by the mass-scaling relation: 

MBH = c|P-P*r, (7) 

where 7 is independent upon the particular choice of the initial data. The nomenclature "type- 
II" comes from the analogous type-II phase transitions in statistical mechanics, which are 
characterized by scale invariance of the thermodynamical quantities [[T3l . 

These phenomena are usually interpreted in terms of attractors in an infinite-dimensional 
phase space, but we will here present a qualitative picture which can be useful to fix the ideas 
(see also the review in flJl ). A more rigorous study employing the renormalization group 
formalism can be found instead in ■ 

Let us consider general relativity as an infinite dimensional dynamical system in an 
abstract phase space in which extra gauge freedoms have been eliminated so that each point, 
Z, can be thought as an initial data-set for the Einstein equations and the associated time 
development as a line in this space: t i-> Z{t). We suppose to have chosen a slicing adapted 
to the self-similarity of the critical solution so that it appears as a fixed point, Z*, in the 
CSS case or a closed orbit for the DSS case (see |fT3l for a more in-depth discussion of the 
consequences of these assumptions). 

In the case of CSS solutions, the main features of this phase space are the presence of 
two attractive sets: A and B representing regular solutions without singularities and black- 
hole solutions. Their basins of attractions are separated by a manifold, C, called critical 
manifold on which there is an attractor of codimension one: the critical solution, Z*; this is 
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shown schematically in figure[T] Any generic one-parameter family of initial data can then be 
thought as a 1 -dimensional line intersecting the critical manifold in one point. Initial data with 
P < P* , will develop as regular solutions not containing singularities and will therefore fall 
in the basin of attraction containing the so called subcritical solutions (cf. set A in figure [T]i. 
Conversely, solutions with P > P* will undergo gravitational collapse with the formation 
of a black hole, thus falling in the basin of attraction containing the so called supercritical 
solutions (cf. set B in figure[T]). 

The key point here is that the critical solution is attractive on the critical manifold. Stated 
differently, nearly-critical solutions will experience "funneling" effects as all but one mode 
converge towards Z*. If P « P*, then the unstable mode, i.e. the mode "perpendicular" to 
C, will be small until later in the evolution, thus allowing for the observation of nearly-critical 
solutions. In this case, all but one mode of the solution are "washed out" by the interaction 
with the critical solution, thus explaining both the universality of the solution and the mass- 
scaling relation. 

2.2.2. Type-I critical phenomena. Type-I critical phenomena are the ones in which the black 
hole formation turns on at finite mass and the critical solution presents a non-selfsimilar 
stationary or periodic solution configuration. The scaling quantity here is the lifetime of the 
metastable solution 

ip = -i In |P-P*|+ const, (8) 

A 

where A does not depend on the initial data. This scaling can be justified using simple 

arguments similar to the ones presented in ifTSl for the mass scaling in the type-II case. 

3. Numerical setup 

In what follows we briefly describe the numerical setup used in the simulations and the 
procedure followed in the construction of the initial data. In essence, we use the Whisky2D 
code described in detail in ||l9l and based on the 3-dimensional code Whisky Ii20ti22l . to 
solve numerically and in 2 spatial dimensions the full set of Einstein equations 

Gf^u ~ StirTfj^i, , (9) 

where G^j^ is the Einstein tensor and T^i, is the stress-energy tensor. More specifically, we 
evolve a conformal-traceless "3+1" formulation of the Einstein equations as presented in 1231 . 
in which the spacetime is decomposed into 3D spacelike slices, described by a metric 7^, 
its embedding in the full spacetime, specified by the extrinsic curvature Kij, and the gauge 
functions a (lapse) and /3* (shift), which specify a coordinate frame. Axisymmetry is imposed 
using the "cartoon" technique ll24l and the equation are solved using finite differencing of 
order three. The chosen slicing condition is the popular "1 + log" while the chosen spatial- 
gauge is the Gamma-freezing one. The field equations for the three-metric 7^ and the 
second fundamental form Kij are coupled with the equations of motion of general relativistic 
hydrodynamics 

V^{pu^) = Q, V^T^^^O, (10) 

where p is the (rest) baryonic mass density, u^ is the four- velocity of the fluid and T^'' is the 
stress-energy tensor of a perfect fluid 

P", = pi?u^u,+p(5^. (11) 
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Table 1. Properties of some of the representative models considered and shown either in 
figure[2]or in figures [6] and [T] More specifically, A^i and Si are the extremes of the range of 
central densities considered. Pi is a largely subcritical model which expands to models P2- 
P4 as the resolution is increased, while Qi and i?i represent the closest super and subcritical 
approximation of the critical solution, respectively. 

Point pc K A/adm Afi, subcritical supercritical 



A^i 


0.00395000 


71.77 


1.3879 


1.5194 


V 


Pi 


0.00459316 


71.39 


1.3832 


1.5194 


^/ 


P2 


0.00341517 


72.23 


1.3754 


1.5077 


^/ 


Pa 


0.00378525 


71.58 


1.3788 


1.5134 


^/ 


Pa 


0.00387685 


71.61 


1.3809 


1.5161 


V 


Qi 


0.00459322 


71.39 


1.3832 


1.5194 




fli 


0.00459322 


71.39 


1.3832 


1.5194 


^/ 


Si 


0.00508840 


71.95 


1.3842 


1.5194 





V 
^ 



Here, H = 1 + e + p/p is the specific enthalpy, p is the pressure, S^^ is the Kronecker 
delta and e is the specific internal energy so that e = p{l + e) is the energy density in 
the rest-frame of the fluid. These equations are closed using an ideal-gas equation of state 
p = (r " l)pe, with adiabatic exponent F = 2. The solution of relativistic hydrodynamics 
equations is obtained via a conservative formulation of ( fTOb as discussed in |fT9l and the use of 
high-resolution shock-capturing (HRSC) schemes with a piecewise parabolic method (PPM) 
for the reconstruction of the primitive variables. The time-stepping is done with a third-order 
total- variation diminishing Runge-Kutta algorithm. Finally, the spatial discretization is done 
on a uniform grid having resolution of either h = 0.1 (medium resolution) or /i = 0.08 (high 
resolution). The outer boundary of the computational domain is set at i? = 15 and we have 
verified that the proximity of the outer boundary does not influence significantly the critical 
solution. 

The equilibrium configuration curves in the (p, A/adm) plane and the perturbative 
oscillations frequencies quoted in the text have been computed using two codes kindly 
provided to us by S'i. Yoshida ll25l and C. Chirenti 



3.1. Initial Data 

The initial data consists of a family of spherical stars having fixed baryonic mass 

A/6 = 1.5194 EE Mb, (12) 

constructed using a polytropic equation of state p = Kp^, with F = 2. Each model is 
computed by fixing its central rest-mass density, pc, while the value of K is fixed after 
imposing the condition (fT2] i. The reason for this choice is that we want to guarantee that 
all the models considered have, at least initially, the same baryonic mass to the precision in 
expression (fTZt . Solutions with different baryonic mass, in fact, are effectively in different 
phase spaces and thus not useful when looking at a critical behaviour. Of course different 
models will also be sUghtly different because the perturbations will slightly alter their mass- 
energy or because although Mb is conserved to high precision by employing a conservative 
formulation of the equations, it is nevertheless not conserved to machine precision. All of 
these latter errors, however, are entirely resolution dependent and can, therefore, be singled 
out by considering simulations at different resolutions. 

These initial models have been evolved under the sole effects of the perturbations induced 
by the truncation error Besides depending on resolution (and converging away), the amplitude 
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Figure 2. Position of some of the most important models in the {pc, Madm) plane, where 
the solid (black) line refers to a sequence with K = 71.39, while the dashed (blue) hne refers 
to a sequence of models having baryonic.mass Mi, = 1.5194 = M^. The points A^i and 
Si are the extremes of the range of central densities considered [cf. eq. 1131 1, Pi is a largely 
subcritical model, while Qi and i?i represent the closest super and subcritical approximation 
of the critical solution, respectively. The inset shows a magnification of the region near the 
critical solution; the properties of the model are reported in table[T] 



of these perturbations is difficult to measure as it depends on a number of different sources 
of error, such as the interpolation error of the one-dimensional initial data on the three- 
dimensional Cartesian grid, or the treatment of low density "atmosphere" regions, which 
are not measurable directly. However, an indirect measure can be obtained by looking 
at a short evolution of a stable spherical star which, in absence of any numerical error, 
would not exhibit any dynamics but which, in practice, oscillates under the effects of these 
perturbations ifTl [191 - 121112714291 . The amplitude of the observed oscillations can be then 
interpreted as an indirect measure of the numerical perturbation. In particular, we can consider 
the value of the average velocity in the radial direction during the first iterations as an estimate 
of the amplitude of an equivalent velocity perturbation. In this case, for a spherical star with 
Pc = 0.00128 and K = 100 evolved for 100 timesteps on a h = 0.1 grid, we measure 
an average velocity, v"^ ~ 1.1 x 10^^. Further insight can also be gained by the average 
of the momentum constraint violation in the radial direction and the Hamiltonian constraint 
violation, which we measure to be ~ 2.3 x 10^^ and ~ 6.1 x 10^^, respectively. 

The determination of the critical value of the central density p* is obtained rather 
straightforwardly via a bisection-like strategy within the initial interval 

0.00395 <pc< 0.0050884 , (13) 

where the extrema correspond to a stable oscillating star or to one collapsing promptly to a 
black hole, respectively. 

The main properties of the initial data are collected in table[T]and summarized in figured 
which reports the position of some of the most important models discussed in this paper in 
the {pc, Madm) plane. More specifically, A'^i and 5*1 are the extremes of the range of central 
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densities considered [cf. eq. (fTTtl. Pi is a largely subcritical model which expands to models 
P2-P4 as the resolution is increased {cf. figure|6]l, while Qi and i?i represent the closest super 
and subcritical approximation of the critical solution, respectively. Note that Ri and Qi differ 
only by the 4.6 x 10~® % in the central density and thus they appear identical in the figure. 
Note also that Pi, Qi and Ri are all on the unstable branch of the models of equilibrium and 
are therefore linerarly unstable. 

As a final remark we note that although the use of an axisymmetric system of equations 
is not strictly necessary for the spherically-symmetric initial data considered here, their 
numerical solutions in 2 spatial dimensions via the Whisky 2D code has been useful in view 
of the connections between the critical behaviour discussed here and the one presented in the 
companion paper Q, where the head-on collision of equal-mass neutron stars is considered. 
The possibility of using the same numerical infrastructure and comparable truncation errors 
has been in fact very important in determining the connections between the two critical 
behaviours. 

4. Results 

In what follows we discuss the nonlinear dynamics of the spherical stars as these evolve away 
from their initial state on the unstable branch and exhibit a critical behaviour. 

4.1. Critical solution 

We first consider the evolution of models in the window ( [13] ) under the sole effect of the 
numerically-induced perturbations. Some of these models, namely the supercritical ones, 
collapse to black hole, while others, namely the subcritical ones, undergo a sudden expansion 
followed by a relaxation towards the corresponding model on the stable branch of the spherical 
star solutions. This is clearly shown in figure[3] which reports the evolution of the central rest- 
mass density and where different lines refer to different initial data in the interval 

0.0045931640625 < Pc < 0.00459371875 . (14) 

By looking at left panel figure [5] it is quite apparent how the survival time of the 
metastable solution increases as the initial models approach the critical threshold and both 
the subcritical and the supercritical solutions overlap for a long part of the evolution, before 
departing exponentially. It is also worth remarking that the linear stability analyses of theses 
models indicates that they are linearly unstable with a characteristic collapse time {i.e. the 
inverse of the imaginary part of the complex eigenfrequency of the fundamental mode) 
T ~ 440. Yet, as shown in figure [3j the metastable models survive for much longer times and 
for almost r ~ 850 for the models closest to the critical threshold. 

A similar behaviour in the evolution of the central rest-mass density has been observed 
also in the simulations reported in |[T|, although those refer to magnetized and rotating stellar 
models and thus, being them result of three-dimensional simulations, are restricted to a much 
smaller interval of significant figures. In addition, and as mentioned in the Introduction, 
evidence for a type-I critical behaviour for the evolution of the central rest-mass density has 
been shown also in the head-on collision of two equal-mass spherical stars ifTOl and will be 
further discussed in the companion paper El. 

As the secular evolution in the central density is a well-known "feature" of the 
numerical solution of relativistic multidimensional stellar models and has been observed 
in codes implementing very different numerical methods and formulations of the Einstein 
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Figure 3. Left panel: evolution of the central rest-mass density near the critical threshold with 
different lines referring to different initial models. Right panel: the same as in the left panel 
but coiTected for the secular evolution given by eq. (15). 



equations |fn [T9]j2ni27H29l . we have isolated this secular behaviour by computing a least- 
square fit of the common part of the evolution in order to isolate the true dynamics from the 
low-frequencies numerical components. More specifically, we have modeled the evolution of 
the central rest-mass density of the metastable equilibrium via the Ansatz 

<Pit) = Po + Pit + P2 cos{2iihit + (y3i) + Pa cos(27r/i2i + if 2) , (15) 

where po ^ P2 are just coefficients in the interpolation and do not have a particular physical 
meaning. On the other hand, the frequencies hi and /12 are chosen as the two smallest 
frequencies appearing in the Fourier spectrum of the central density during the metastable 
phase {cf. figure|4]and see also discussion below on the spectral power density of the putative 
critical solution). The residuals after the fit are shown in the right panel of figure |3] and help 
considerably in appreciating the dynamics of the unstable models near the critical value. 

Using a large set of simulations with resolution of h = 0.1 and a straightforward 
bisection strategy we have located the critical threshold to black-hole formation at a central 
density 

p* = 0.004593224802 ± 2.1 x 10"^^ . (16) 

Clearly, we expect this value to depend on the initial perturbation and thus on the resolution 
used, as well as on the numerical method employed. On the other hand, we also expect 
that the associated solution and the critical exponent to be "universal", in the sense that they 
should not depend depend sensitively on the perturbation or on the particular family of initial 
data as far as this family is characterized by a single parameter and thus intersects the critical 
manifold C in a single point which is near enough to this solution. In this case, in fact, the 
associated critical solution is supposed to be at least locally attractive on a sub-manifold of 
the phase space of codimension one. 

To validate that the behaviour discussed so far and shown in figure [3] does represent a 
type-I critical behaviour we compute the survival time of the metastable solution r, i.e. the 
"escape time", and study how this varies as the critical solution is approached. We recall that 
we expect that the escape time near the critical for a type-1 critical phenomena should behave 
as 



T = -^\ll\pc- Pa\ 



const , 



(17) 
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Figure 4. Power spectral density of the evolution of the central rest-mass density for the model 
closest to the putative critical solution (i.e. with pc = 0.0045932248034) when the secular 
drift part 1151 has been removed from the data. The eigenfrequencies associated with the 
corresponding spherical star model are also shown as vertical lines. 



and such expected solution is indeed shown as a dashed line in figure |5] Also shown with 
squares and triangles are the computed escape times for different initial data and different 
resolutions (blue squares for h ~ Q.l and red triangles for h ~ 0.08). The latter are calculated 
in terms of the time t^ at which the relative difference between the observed central baryonic 
density and the best approximation of the critical solution <f>{t) ( fTSl ) becomes larger than e. 
We find that, for a large enough e, such that 1 ^ e > e* > 0, these times depend only weakly 
on e and thus give a good measure of the departure time from the critical solution. A value 
of e = 0.5% provides a sufficiently accurate measure and this is the one employed for the 
data points shown in figure |5] We finally estimate the critical exponent A by making a linear 
least-square regression of the data points of sub- and supercritical solutions and then by taking 
the average of the two values. Using the medium-resolution h = 0.1 simulations we therefore 
obtain for the critical exponent 

A = 0.02149665, (18) 

with a coefficient of determination R^ relative to the linear regression ( fTTI i and computed 
on the full dataset containing both sub and supercritical solutions, of 0.960517. The critical 
exponent ( fTSl l is found also in the case of the h ~ 0.08 simulations, although in this case 
the scattering is somewhat larger and the data agrees within 7%. We note that these high- 
resolution simulations are computationally very expensive and this is why we have restricted 
them to a smaller set of initial data. Clearly, the match between the computed escape times 
and the one expected from the critical behaviour is very good over the 6 orders of magnitude 
in \pc — p*\ spanned by our data-set and thus provide convincing evidence that indeed critical 
behaviour can be found in the dynamics of linearly unstable spherical stars. 

As a final remark we note that while the evidence for a critical behaviour is clear, much 
less clear is the physics of the critical solution which is, after all, a perturbed spherical star 
Recent studies of nonlinear perturbations of relativistic spherical stars have shown that linearly 
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Figure 5. Escape time t as a function of In | p — p* | for subcritical (left panel) and supercritical 
solutions (right panel), respectively. The blue squares corresponds to the results obtained with 
the h = 0.1 resolution, while the red triangles to the results obtained with the h = 0.08 
resolution. The dashed lines represent the fit obtained using I17t with the A obtained from the 
h = O.l solutions. 



unstable stars can be stabilized via nonlinear couplings among higher-order modes OOll . It 
is possible that such a nonlinear coupling is present also here and we conjecture therefore 
that the stability of the metastable solution is due to mode coupling of the first overtones of 
the fundamental mode. Support to this conjecture comes from the power spectral density in 
figure m which shows that, apart from the F-mode which is obviously missing as it has only 
imaginary eigenfrequency, the spectrum of the metastable solution is essentially identical to 
the one of an excited spherical star with {pc, K) ~ {p* ,K*) and M{, = Mf,. Interestingly, 
most of the energy is in the first overtone, HI, even though the numerical perturbation can 
be thought as "white noise" exciting all the modes of the star with almost equal energy. 
The behaviour discussed above persists also when considering models with higher spatial 
resolutions. 



4.2. Subcritical solutions 

While the final fate of supercritical solutions is clearly that of leading to a collapse and 
to the formation of a black hole, the one of subcritical solutions deserves a more detailed 
explanation. As one would expect, given that the initial data represent linearly unstable 
stars, the subcritical solutions show a first expansion as the star migrates to the stable branch 
of the equilibrium configurations, which is then followed by a slow relaxation where the 
central rest-mass density exhibits strong oscillations around smaller and smaller values, that 
would eventually reach in the continuum limit, the value corresponding to the model on the 
stable branch having the same gravitational mass of the initial one. In practice, however, the 
migration to the stable branch is accompanied small losses both in the gravitational mass and 
in the rest-mass which, although smaller than ~ 0.7%, need to be taken properly into account. 
More specifically, we have analyzed in detail the evolution of the largely subcritical 
model Pi, {cf. table [T]), which is an unstable spherical star with an F-mode whose imaginary 
part of the eigenfrequency is i^i = 0.461 kHz. We evolve therefore evolved such a model 
it with three different spatial resolutions of h = 0.1, h = 0.09 and h ~ 0.08, and studied 
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its migration to the stable branch. The asymptotic state of the solution and in particular to 
the final central rest-mass density pj is estimated by modelling the time evolution of the 
oscillating star on the stable branch with a simple Ansatz of the type p{t) ^ pj: + p^/t and 
by performing a nonlinear least square fit on an appropriate window including the final part 
of the dynamics. For any given resolution we have then computed the total baryonic-mass 
losses due to the numerical dissipation AMb = Mb — Mbj, and determined the polytropic 
coefficient Kf yielding a spherical stellar model with central rest-mass density pf and baryon 
mass Mb J. Clearly, for such a model it is then also possible to compute the gravitational mass 
and thus track the migration on a {pc , Mj^^-^^ ) plane. 
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Figure 6. Dynamics of the migration on a {pc , J^^adm ) plan^- The hneai'ly unstable and 
largely subcritical model Pi migrates to a new solution on the stable branch of equiUbrium 
configurations. Indicated with P2 — Pi are the new asymptotic states for resolutions h = 
0.1 — 0.08, respectively. Indicated with a thick solid line is the sequence of initial models 
having the same polytropic index of Pi, while indicated with dotted and dashed lines are 
sequences of models having the same rest-mass as the asymptotic models P2 — P4. Finally, 
shown as Pi is the asymptotic state of Pi in the continuum limit; note that even for the coarse- 
resolution case the changes in baryonic and rest-mass are only of ~ 0.7%. 



The overall results of these migrations are shown in figure|6l where we report the stellar 
configurations on Alb — const, curves. The minimum of each curve corresponds to the 
maximum in the usual {p, Madm), K = const., plots and separates the stable and unstable 
branches of solutions. When a resolution of h = 0.1 is used the model Pi migrates to the 
new asymptotic model P2, while it will migrate to models P3 and P4 as higher resolutions of 
h = 0.09 and h = 0.08 are used, respectively. Note that already with the coarsest resolution 
of h = 0.1 the losses in gravitational masses are ~ 0.65% and that these decrease to ~ 0.16% 
when a resolution of h ^ 0.08 is used. Finally, indicated with Pi is the expected asymptotic 
model when the numerical losses are extrapolated to the continuum Umiyl; clearly, in the limit 
h -^ 0, the migration of model Pi takes place to a new state having the same gravitational 

f Note that we do not mark this point with a symbol as it does not correspond to a numerically computed value, as 
it instead for P2, P3 and P4 
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and baryonic mass as the initial one. 

4.3. Perturbation of nearly-critical solutions 



As discussed in Sect. 14.11 the central rest-mass density of the linearly unstable models can be 
used as a critical parameter for the gravitational collapse of a linearly unstable spherical star, 
in contrast to what has been observed for example by Novak in Q or by Noble in ID . We 
believe this is due to the very different set of initial data selected here and in IJUS). Indeed, the 
reason why this behaviour has not been observed in many previous studies is that we consider 
initial stellar models that are already linearly unstable, in contrast with what done in lllll], 
where the initial models are instead linearly stable and then subject to a perturbation (either by 
introducing a radial velocity Ii2j|9], or by considering employing the interaction with a scalar 
field IS)). For our set of initial data, therefore, the critical solution is essentially a spherical 
star with an unstable F-mode, and any finite perturbation exciting this mode will change the 
solution in a dramatic way (A discussion of this change within a phase-space description will 
be made later on when presenting figure [8]1§- 

To confirm this hypothesis, we follow 121 and IH, and construct a new family of spherical 
initial data obtained by perturbing the slightly supercritical model Qi {cf. table [T]) via the 
addition of a radial velocity perturbation in the form of the 3-velocity component 

i;'-(x) = ^(3x-x3), 2;=-^, (19) 

Z Hi, 

where U is the amplitude of the perturbation at the surface of the star, i?* and can be 
either positive (outgoing radial velocity) or negative (ingoing radial velocity). Because the 
perturbation ( fT9l ) matches the eigenfunction of an idealized F-mode perturbation, it should 
excite the only unstable mode of the critical solution. 

Performing simulations for different values of U and a resolution h = 0.1 we 
find, not surprisingly, that for negative values of U the perturbed models of Qi collapse 
to a black hole. Furthermore, because in this case the radial velocity accelerates the 
development of the unstable mode, the larger the values of U the shorter the time to collapse, 
i.e. T ^ — ci log(C/) + C2, where ci and C2 are positive constant coefficients. On the other 
hand, for positive values of U, the perturbed models of Qi, which we recall is supercritical 
for U = 0, becomes subcritical and shows the same qualitative behaviour as that of model 
i?i. Hence, a suitably perturbed supercritical model can behave as a subcritical one. 

The dynamics of these perturbed, nearly-critical models is shown in figure [T] where 
the solid (black) line represents the supercritical solution Qi, while the dotted (blue) line 
represents the subcritical solution i?i. The dashed lines show again the evolution of Qi, 
but when subject to a positive (red short-dashed line) or negative (green long-dashed line) 
velocity perturbation. The dynamics shown in figure |7] underlines an important characteristic 
of critical phenomena: the precise value of the critical parameter at the intersection between 
the one-parameter family of solutions and the critical manifold depends on the family itself. 
In particular this means that if we fix a value of the perturbation amplitude, [/ / 0, we have 
to expect to find the critical solution at a value of p*{U) different from the one quoted in ( fT6b 
which is attained in the case C/ = 0. For this reason the application of a non infinitesimal 
perturbation to a nearly-critical solution results in a dramatic change in the dynamics of the 
system. 

§ With "perturbation" we are here referring to a globally coherent, resolution independent perturbation such as the 
one given in eq. {\9\ . This has to be contrasted with the random, truncation-error induced and resolution-dependent 
perturbations we have considered in Sect. 14. II 
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Figure 7. Perturbation of nearly-ciitical solutions. The solid (black) line represents the 
supercritical solution Qi, while the dotted (blue) line represents the subcritical solution Ri. 
The dashed hues show again the evolution of Qi, but when subject to a negative (red short- 
dashed Une) or positive (green long-dashed line) velocity perturbation. Clearly, in the latter 
case the supercritical solution Qi becomes subcritical and shows the same behaviour as the 
solution Ri . 



The phase-space representation of this concept is summarized in figure[8] where we show 
two one-parameter famines of perturbed TOV initial data, whose critical parameter, p^, is the 
central rest-mass density. The perturbation is given by the composition of truncation errors 
and of a radial velocity perturbation U in the form ( fT9l ), where U = or U ~ Uq > 0. As 
these families represent different initial configurations, they will intersect the critical manifold 
C at two different points, with correspondingly different values of the critical parameter 
{0, /o*(0)} and {Uo, ptiUo)} (these points are marked as filled circleslU|. In particular, when 
U runs between and Uq, the set of critical configurations {U,p*{U)} will represent a curve 
on the critical manifold C and this is shown with a violet solid line in figure [8] Considering 
now a configuration near {0, p* (0) } and applying to it a velocity perturbation in the form ( fT9l ) 
with U = Uq, will produce a new configuration {Uq, p*(0)} which is not necessarily on the 
critical manifold (this is marked with a filled square). Indeed, the whole family {U, p*(0)}, 
that is the set of configurations with a nonzero initial velocity perturbation but central density 
which is the critical one for the zero-velocity case, are in general expected to be outside the 
critical domain. The family {U, p* (0)} is shown with a black dot-dashed line in figure[8] 

As a final remark we note that another important difference between the work presented 
here and that in |l7]|9l is that the we find evidence of a type-I critical behaviour with a periodic 
solution, in contrast to what found in II7I9II . which is instead of type-II and with DSS solutions. 
We believe the origin of this important difference and of the presence of a periodic solution is 
in our use of an ideal-fluid EOS and hence in the presence of an overall scale in the problem. 

II In our notation, the point {Uq, p'^{Uo)} is the ciitical solution with initial velocity perturbation given by 1191 with 
U = Uo- Similarly, a configuration {Uo, pj:(0)} will be a member of the family with initial velocity perturbation 
C/q , but with a central density which is the critical one for a model with U = 
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Figure 8. Phase-space diagram representation of nearly-critical solutions. In particular, we 
show with red solid curves two one-parameter families of initial data, perturbed by the addition 
of a radial velocity profile in the form U9) with U = or U = Uo > 0. The locus of the 
critical points, {U, Pc{U)}, is shown with a violet solid line, while the family of initial-data 
{U, p* (0)} is shown with a black dot-dashed hne and the point {Uq, p* (0)} is marked with 
a filled square. The latter represents therefore the family of initial data obtained by adding a 
velocity perturbation with amplitude U to the model with central density would when U = 0. 
Also highlighted with filled circles are the critical points for the families with (7 = and 
U = Uo, i.e. {0,p*(0)} and {I/q, P*(C/o)}. 



Conversely, the spherical stars considered in the above mentioned works were evolved using 
either an ultrarelativistic EOS |5] (which, as commented in the Introduction, are intrinsically 
scale-free) or with very strong perturbations ||7][9l, thus in a regime of the EOS which is 
approximatively ultrarelativistic |fT3l . 



5. Conclusions 

In general, critical phenomena in gravitational collapse are of great interest because they play 
a central role in phase transition of families of solutions in general relativity. In more specific 
context of the dynamics of NSs, type-I critical phenomena have seen a renewed interested 
when it was shown that a critical behaviour of this type is produced in the in head-on collision 
of NSs ifTOl or in the dynamics of rotating magnetized stars fT\. With the goal of studying in 
more detail the occurrence of type-I critical collapse in NSs, we have therefore employed the 
2D general relativistic code Whisky 2D to study a large set of spherical stellar models having 
a constant baryon mass. Differently from what done before by other authors, e.g. EHH, we 
have considered stellar models that are on the "right" branch of the models of equilibrium and 
thus linearly unstable. 

Using a simple ideal-fluid EOS and very small perturbations which are entirely induced 
by the truncation error, we have found that our family of initial data exhibits a clear type-I 
critical behaviour at at a threshold central rests-mass density of p* ~ 0.004593224802 ± 
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2.1 X 10^^^ and with a critical exponent A = 0.02149665. These results thus confirm the 
conclusions reached by Liebling et al. [Tj but also provide a more quantitative determination 
of the threshold and of the nature of the critical scaling. Exploiting in fact the relative 
simplicity of our system, we were able carry out a more in-depth study providing solid 
evidences of the criticality of this phenomenon and also to give a simple interpretation of the 
putative critical solution as a spherical solution with the unstable mode being the fundamental 
F-mode. As a result, we have shown that for any choice of the polytropic constant, the 
critical solution distinguishes the set of subcritical models migrating to the stable branch of 
the models of equilibrium from the set of supercritical models collapsing to a black hole. 

Furthermore, we have studied how the dynamics changes when the numerically 
perturbation is replaced by a finite-size, resolution-independent velocity perturbation and 
show that in such cases a nearly-critical solution can be changed into either a sub or 
supercritical. Finally, the work presented here here is of direct help in understanding why 
the critical behaviour shown in the head-on collision of two neutron stars is indeed of type-I 
and why it can be explained simply in terms of the creation of a metastable stellar model on 
the unstable branch of equilibrium solutions . 
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